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Abstract 

We consider low-energy nucleons at next-to-next-to-leading order in lattice chiral effective 
field theory. Three-body interactions first appear at this order, and we discuss several 
methods for determining three-body interaction coefficients on the lattice. We compute 
the energy of the triton and low-energy neutron-deuteron scattering phase shifts in the spin- 
doublet and spin-quartet channels using Liischer's finite volume method. In the four-nucleon 
system we calculate the energy of the a particle using auxiliary fields and projection Monte 
Carlo. 
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I. INTRODUCTION 



We study low-energy nucleons on the lattice at next-to-next-to-leading order in chiral 
effective field theory. In Weinberg's scheme counting orders in effective field theory 

is equivalent to dimensional analysis for irreducible diagrams. The expansion parameter is 
Q/ A, where Q is the low momentum scale associated with external nucleon momenta or the 
pion mass, and A is the high momentum scale at which the effective theory breaks down. 
Terms at next-to-next-to-leading order are of size Q^/A^, and three-nucleon interactions 
first contribute at this order. In this work we consider three-nucleon forces on the lattice 
for systems with three and four nucleons. Our analysis continues a series of recent papers 
on lattice chiral effective field theory for few- and many-nucleon systems. Previous studies 
lave considered dilute neutron matter and light nuclei using interactions at leading order 

Our discussion is organized into three parts. The first part begins with an overview 
of the effective potential for nucleons up to next-to-next-to-leading order in chiral effective 



field theory. Reviews of chiral effective field theory can be found in Ref. 



Id, 



nj. We 



discuss some simplifications that can be made at low cutoff momentum and present lattice 



operators for each interaction. Nucleon-nuc 



determined using a spherical wall method 



eon phase shifts and the S-D mixing angle are 
12| and used to set unknown operator coefficients. 



In the second part we calculate the low-energy spectrum for three nucleons. We compute 
the triton energy and determine neutron-deuteron phase shifts using Liischer's finite volume 
method. These are used to constrain the two unknown three-body operator coefficients. 
In the third and final part we rewrite the lattice action in terms of auxiliary fields and use 
projection Monte Carlo to calculate the energy of the a particle. This leads to a discussion 
of alternative methods for fixing three-body operator coefficients. We summarize our results 
and discuss possible extensions in future work. 

II. CHIRAL EFFECTIVE FIELD THEORY 
A. Effective potential for two nucleons 

In the following q denotes the t-channel momentum transfer for nucleon-nucleon scattering 
while k is the w-channel exchanged momentum transfer. We assume exact isospin symmetry 



and neglect electromagnetic interactions. At leading order (LO) in the Weinberg scheme the 
two-nucleon effective potential consists of two independent contact terms and instantaneous 
one-pion exchange (OPEP), 

V^o = V^^^ + yOP^P. (1) 

The scattering between nucleons consists of contributions from direct and exchange dia- 
grams. Nevertheless for bookkeeping purposes we label the interactions according to the 
tree-level scattering amplitude for distinguishable nucleons. For two-nucleon interactions 
we label one nucleon as type A and the other nucleon as type B. In this notation the 
amplitude for V^^^-* is 

^[y(o)] =Cs + CT{aA-aB). (2) 

and the amplitude for \/Opep 

The vector arrow in a signifies the three-vector index for spin. The boldface for r signifies 
the three-vector index for isospin. We take for our physical constants m = 938.92 MeV 
as the nucleon mass, = 138.08 MeV as the pion mass, = 93 MeV as the pion decay 
constant, and qa = 1-26 as the nucleon axial charge. 

At next-to-leading order (NLO) the two-nucleon effective potential contains seven in- 
dependent contact terms carrying two powers of momentum, corrections to the two LO 
contact terms, and the leading contribution from the instantaneous two-pion exchange po- 
tential (TPEP) 3, 14, 15, 30, 

Fnlo = ^Lo + Al^(°) + 1^(2) + V™P. (4) 

The tree-level amplitudes for the contact interactions are 

A [AV(o)] = AC5 + ACt (c?a ■ ^b) (5) 

and 



+ Ce {a A ■q){aB-q) + Cj [a a ■ (<Tb ■ Pj ■ (6) 
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The amplitude for the NLO two-pion exchange potential is 



3847r2/. 



^ml (5,1 - 4,1 - 1) + (23,1 - I0g\ - l) + 



647r2/, 



where 



Tf \ ^ I A 2,2 1 V4m2 + g2 + g 



(7) 



(8) 



At next-to-next-to-leading order (NNLO) there are no additional two-nucleon contact 
interactions, but the two-pion exchange potential contains a subleading contribution, 

3,3 



A r;/TPEPi 

L^^NNLO J 



167r/; 



- A{q) {2ml + 1^) [2m^ (2ci - C3) - Cgg'] 



A(g) (4m2 + [(g-. 5^^) (g-. _ ^2 . 



where 



. , . 1 g 

AiQ) = — arctan 

2g 2m^ 



(9) 



(10) 



The low-energy constants ci, 03,04 parameterize the coupling of the nucleon to two pions. 
These constants have been determined from fits to low-energy pion-nucleon scattering data 
[20I . and in the following we use the values ci = —0.81 GeV~\ C3 = —4.7 GeV~\ C4 = 



3.4 GeV-i 



2l|. 



B. Three-nucleon interactions 



A number of different phenomenologica 



the literature 



22 



23 



24 



25 



26 



MM 



;hree- nucleon potentials have been introduced in 
□ 

,|30|. Effective field theory provides a systematic 
method for estimating the relative importance of three-body interaction terms. Few-nucleon 
forces in chiral effective field theory beyond two nucleons were first introduced in Ref . 2[ . In 
Ref. 311 it was shown that three-nucleon interactions at NLO cancel and three-body effects 
first appear at NNLO. The NNLO three-nucleon effective potential includes a pure contact 
potential, v}^^l^^., one-pion exchange potential, , and a two-pion exchange potential. 



(3JV) 



TPE ' 



(37V) 



NNLO 



(37V) 



^ ^OPE ^ 



(37V) 
TPE ■ 



(11) 



The corresponding diagrams are shown in Fig. [TJ 



/ \ / 












Tl 


n 



(a) 



(b) 



(c) 



FIG. 1: Three-nucleon forces at NNLO. Diagrams (a), (b), and (c) show the contact potential, 



(3JV) 



, one-pion exchange potential V^^J , and two-pion exchange potential Vr 



{37V) 



TPE • 



Similar to our bookkeeping notation for two-nucleon interactions, we write the tree-level 
amplitude for three-nucleon interactions where the first nucleon is type A, the second nucleon 
type B, and the third type C. We sum over all permutations P{A, B, C) of the labels, and 
^4; are defined as the differences between final and initial momenta for the respective 



nucleons. The amplitudes for V^l^2.ct ^-nd Vqpe' are 



.{37V) 



32 



33| 



A 



{37V) 



contact 



A 



Vr 



(37V) 



OPE 



P(A,B,C) 
QA ■ ^A 



9A 



{qA ■ <?i?) {ta ■ Tb) ■ 



(12) 



(13) 



The coefficients E and D are both cutoff dependent. The coefficient E determines the 
short distance interactions between three nucleons, while D determines the pion coupling 
to two nucleons. Following the notation introduced in Ref. 
parameters ce and cd such that 

ce CD 



33| . we define dimensionless 



E 



D 



700 MeV. 

For convenience we separately label three parts of the two-pion exchange potential. 



(14) 



where A^^ ~ nip. We take A^^ 



(37V) 



TPE 



. .{37V) . .(37V) . .(37V) 
*^TPE1 "T '^TPE2 "T *^TPE3- 



TPE2 



(15) 



The corresponding amplitudes are 



A 



(37V) 



TPEl 



C3 / ^ \ 
/2 ^2/. i ^ 



U ■ 



^a) {qb ■ ^b] 



P{A,B,C) 



{q\ + ml) {q% + ml] 



{Qa ■ Qb) {ta ■ Tb) 



(16) 



A 



(37V) 



TPE2 



J 7r 



2/. 



>p {qa ■ ^a) {qb ■ ^b) 



ta ■ Tb) 



(17) 



(37V) 



TPE3 



C4 



2/2 V2/. 



9A 



E 



P(A,B,C) 



(gA • ^a) (qb ■ ^b) 

{q\ + ml) {ql + ml] 



[{qA X Qb) ■ ^c] [(ta X Tb) ■ Tc] 



(18) 



C. Simplified form at low cutoff momentum 



A number of alternatives to Weinberg's power counting have recently been discussed in 
the literature |3J, |35|, l36|, l37|. However at low cutoff momentum the advantages of these 
alternative schemes are numerically small [38]. In this study we use spatial lattice spacing 
a = (100 MeV)~^, corresponding with cutoff momentum A = 314 MeV ~ 2.3m7r. Our choice 
of low cutoff scale avoids numerical problems in Monte Carlo simulations due to spurious 
deeply-bound states and large sign/phase oscillations. 

In the following lattice calculations we use Weinberg's power counting with some addi- 
tional simplifications made possible by the low cutoff momentum. For nearly all |g| < A we 
can expand the NLO two-pion exchange potential in powers of q'^/{Aml), 



A[V, 



TPEPl 
NLO J 



ta ■ Tb 



3847r2/.^ 



4ml {8g\ - - 1) + (3%1 - I7g\ - 2) + mlO 



q 

im'i 



[{q- ^a) (g ■ ^b) - {^A ■ ^b)] 



(19) 



This expansion fails to converge only for values of q near the cutoff scale A 2.3m^, where 
the effective theory is already problematic due to large cutoff effects. From a practical 
viewpoint there is no advantage in retaining the full non-local structure of V^lq"^ at this 
lattice spacing. Instead we simply use 



NLO 



Vio + A\/W + \/(2) 



(20) 
(21) 



where the terms in Eq. f[T^ with up to two powers of q are absorbed as a redefinition of the 
coefficients Ay(o) and V'^^l 
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Similarly we can expand the NNLO two-pion exchange potential, 



A [T/TPEPl 



■ 167r/4 



(2ci - cs) + (^c, - He3) + mlO ((4)') 



[(<? ■ o'a) (g ■ o'b) - {^A ■ Ob 



(22) 



327r/4 

The terms with two powers of q were already included at NLO, and so there are no addi- 
tional terms in the two-nucleon potential at NNLO. In our low cutoff scheme the only new 
contributions at NNLO are due to three-nucleon interactions, 



^NLO — ^LO + ^NLO- 



(23) 



III. LATTICE INTERACTIONS AT LO AND NLO 



A. Transfer matrix at LO 



In our Euclidean-time lattice formalism the transfer matrix operator is the normal-ordered 
exponential of the lattice Hamiltonian, : exp(— At) : , where At equals one temporal lattice 
spacing, a^. At leading order we use the LO2 transfer matrix with Gaussian-smeared 

Since we consider only one action, we drop the "2" subscript on LO2. 
The transfer matrix operator is 



Ml 



LO 



: : exp <j -i/freea* ~ ^^fi^ 



Cp^'''^iq)p'^'''^i-q) + CpY.pi^\q)pi^\-q) 



+ 



XI 5Z '^s^S2{ni - n2)ps,'}ini)ps^}{ 



n2 



(24) 



81,82,1 ni,n2 

The momentum-dependent coefficient function f{q) is given by 



f{q) = fo ^ exp 



where 



exp 



-"Ed 



-''Ed 



cos qi ) 



cos qi^ 



(25) 



(26) 



Our lattice notation is defined in the Appendix. The densities p'^^'"' and pf'"" and spin- 
dependent one-pion exchange potential Gs^Sj also defined in the Appendix. We use 



non -perturb ative 
"improved" LO 



perturb ative 
expansion in QIK 




FIG. 2: The "improved" LO action is iterated non-perturbatively while the remaining higher-order 
interactions are treated using perturbation theory. 



the value h = 0.6, which gives approximately the correct effective range for the two S-wave 
channels when C and Cp are tuned to the physical scattering lengths. C is the coefficient of 
the Wigner SU(4)-invariant contact interaction 39|], and Cj2 is the coefficient of the isospin- 
dependent contact interaction. In terms of coefficients for the isospin-singlet and triplet 
channels, 

C = (3C^=^ + C^=°) /4, (27) 

Cp = {C'=^ - C^=o) /4. (28) 

This "improved" leading-order action is treated non-perturbatively while higher-order in- 
teractions are included as a perturbative expansion in powers of Q/A. This is sketched in 
Fig. El 

In pionless effective field theory the three-nucleon contact interaction is included at lead- 



ing order 



4G 



41 



42]. This is needed to stabilize the three-nucleon system in the limit 



of zero-range interactions 



431]. In this study we use chiral effective field theory where the 



interactions have nonzero range. However given our coarse lattice spacing, we may find 
that the three-body contact interaction is numerically large and requires non-perturbative 



treatment. Non-perturbative t 
cussed in the literature 



44 



45 



iree-body contact interactions on the lattice have been dis- 



461 . However for the lattice calculations presented here we 



choose a different approach to address the same problem. 

We use the fact that the three-nucleon interaction depends on both the spatial lattice 
spacing, a, and the temporal lattice spacing, a^. The temporal lattice spacing regulates 
the transfer matrix element when the interaction potential energy exceeds a^^. As a result 
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it can affect the magnitude and sign of the three-body contact interaction [47]. With the 
spatial lattice spacing held fixed, we dial the temporal lattice spacing to a value where the 
three-nucleon interaction is numerically small. This involves a calculation of the spectrum 
of the three-nucleon system. In the following calculations we use the lattice spacings 
a = (100 MeV)^^ and at = (150 MeV)^^ and show that the strength of the three-nucleon 
contact interaction is small enough to be treated using perturbation theory. For these 
lattice spacings we find leading-order coefficients C^^^ = —5.105 x 10~^ MeV~^ and C^^^ = 
—3.507 X 10^^ MeV^^ when tuned to the physical S-wave scattering lengths. 

B. Transfer matrix at NLO 

At next-to-leading order the lattice transfer matrix is 



Mnlo = Mlo - at : [AV + AVj2 + V,2 + Vp^g2 + Vs2,g2 

+Vs^,p,q^ + V(g.5)2 + Vi2^(^g.s)2 + V({=x5).fc] Mlo : . 



(29) 



The corrections to the leading-order contact interactions are 




(30) 



n 




(31) 



and the seven independent contact interactions with two derivatives are 




(32) 



n,l 




(33) 



n,I,l 




(34) 



n,S,l 



V52,/2,q2 = ~Cs2j2^g2 ! ^ pt J {n)^^ f {u) !, 



(35) 



n,S,I,l 




-s)2 = \C{qsy ■^^Aspi'''{n)^As^pi,'''{n) :, 



(36) 



n S S' 
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(37) 



n,I S 



YI=1 



I 

— ( 

2 



n,LS,l' 



Ilf'%n)A,p^J'%n) + </(n)A,p'^^'^(n) 



(38) 



n,l,S,l',I 



The densities, current densities, and symbols A/ and V;^, are defined in the Appendix. 

The Vf'^^^gyf, term is designed to ehminate lattice artifacts in the spin-triplet even-parity 
channels. This is done by projecting onto the isospin-triplet state. 



(39) 



where 

y{iqxS)-k 
Vl'2,{iqxS)-k = 

and 



--C{iqxS)-k '■ ^ £l,S,l' 
n,LS,V 



nf^'^{n)Avpf\n) + </(n)A,p-^'^(n) 



(40) 



-C 



P,{iqxS)-k 



n,l,S,l',I 



^l,S,l' 



</(r!)A,p^;f(n) + nf,}(n)A,pf''^(rI) :, (41) 



{iqxS)-k — ^^(iqxS)-k} 

n — ^ nT^=\- 

'^I^,(iqxS)-k — ^'-^(igx5)-fc- 



(42) 
(43) 



We measure phase shifts on the lattice by imposing a spherical wall boundary on the 
relative separation between two nucleons at some chosen radius. From the properties of the 
spherical standing waves we determine scattering phase shifts and mixing angles 12j. The 
scattering results are nearly identical with the LO2 data at lattice spacings a = (100 MeV)~^ 
and at = (70 MeV)~^ presented in Ref. [5j]- The values for the next-to- leading order 
coefficients are shown in Table [B These values are similar to the NLO2 coefficients in 
Table III of Ref. [St], though there are some differences due to the change in temporal lattice 
spacing. 
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TABLE I: Results for the NLO operator coefficients 



Coefficient 




AC fMeV^l 


4.08 X 10"^ 


ACj2 [MeV^^] 


5.92 X 10^^ 


C 2 fMeV"^l 


-1.31 X 10"^ 


Cr2 2 fMeV""^! 


-3.26 X 10"^° 


C52^,2 [MeV-4] 


-1.53 X 10"^° 


C52,,2,,2 [MeV-4] 


-2.64 X 10-1° 


q,.5)2 [MeV-4] 


-1.92 X 10^1° 


C,2,(,.s)2 [MeV-4] 


9.20 X 10-12 




1.11 X 10-1° 



IV. THREE-NUCLEON INTERACTIONS AT NNLO 



At next-to-next-to-leading order the transfer matrix is 



(SiV) 



.(37V) 



contact 



(37V) 



.(37V) 



'OPE 



'TPEl 



TPE2 



'TPE3 



Ml 



LO 



(44) 



From the constraints of isospin symmetry, spin symmetry, and Fermi statistics, there is only 
one independent three-nucleon contact interaction 33|, |42| . In Eq. (fT2|) we wrote this as a 
'Ta ■ 'Tb interaction over all permutations of the labels A,B,C. For our lattice action we 



choose to write the contact interaction v}^^^^^ as a product of total nucleon densities. 



yi3N) _ i r) 
•^contact g -^contact 



n 



(45) 



The one-pion exchange potential Vqj^ can be written as 



n,S,In',S' 



The three two-pion exchange terms V^p^i, V"r^PE2' ^pes ^'^^ 



(37V) 



TPEl 



X (A5"VrKn^nJA57r;(n,nO) :p^!'?(n')P5"/(^'')p'^''''(^) : 



(46) 



(47) 
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l4S = Wmj||i5:^^ 

n,I n',S'n",S" 



{As'n'j{n, nt)n']T'j{n, rit)) 



X {/\s'<{n\nt)UT,',{n,nt)) : pi'}{n')pifAn")p'''^\n) : 



(48) 



(37V) 



TPE3 - 



(49) 



Definitions for tlie □ symbol and tlie two-point pion correlation functions are given in the 
Appendix. 

In the continuum limit the tree-level scattering amplitudes are 



A 



y{3N) 1 _ r) 

contact -L^ contact ; 



A 



Vr 



(37V) 



OPE 



-D 



OPE 



9a Qa ■ o'A 

•"^ P{A,B,C) 



7a ■ (^B (ta ■ Tb 



A 



(37V) 
TPEl 



D 



9\ 



TPEl 



{qa ■ oa) {qb ■ ^b) 
4/^ P(i;^,c) + "^-^ ^^^^ + 



E 



(^4 • gs) {ta ■ Tb) 



A 



(37V) 



TPE2 



{qa ■ 5a) {qb ■ ^b) 



[Ta ■ Tb) 



(50) 
(51) 

(52) 

(53) 



A 



(37V) 



TPE3 



^TPE3f^ 



{qA ■ ^a) {qb ■ 5b) 
'4/^ p^'x^^c^ (Qa + K) ill + 



[iu X Qb) ■ ac] [{ta X Tb) ■ tc] ■ 



Comparing these with Eq. fll2til8p . we have 

3ce 



D 



TPEl 



-3E 



C3_ 

P'' 
J n 



DtpE2 ■ 



D 



CD 



OPE 



2c, 

J vr 



4/. 4/3A,' 
2/f 



^nr, -L^TPES — :r^- 



(54) 

(55) 
(56) 
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FIG. 3: Plot of CE versus cd when constrained by the triton energy at infinite volume. The solid 
line shows data matched to the physical triton energy, —8.48 MeV. The dotted line shows results 
fitted to the pseudo triton energy, —8.68 MeV, estimated by replacing nn forces with np forces. 



A. Triton energy 

With the lattice transfer matrices Mlo, Mnlo, ^nnlo, we use iterative sparse-matrix 
eigenvector methods to compute the triton energy for cubic periodic lattices. We consider 
cubes with side lengths L = 3, 4, 5, 6, 7, 8 and extract the infinite volume limit using the 



asymptotic parameterization 



c 

-E'triton(-^) ~ -E'triton " -^/^O _ ^g-J?^ 



Lq is a length scale associated with the physical size of the triton wavefunction. For the 
NNLO calculation we fix the coefficient ce as a function of cd by matching the physical triton 
energy at infinite volume, —8.48 MeV. This constraint produces the solid line shown in Fig. 
|3l In the same figure the dotted line shows results obtained by fitting to the "pseudo" 
triton energy, —8.68 MeV. This pseudo energy is an estimate of the triton energy when 
nn interactions are replaced with np interactions [33] • This adjustment takes into account 
the systematic error in our isospin-symmetric calculations with two-nucleon interactions 
matched to np phase shifts. We note the similarity between Fig. [3] and other plots of 



versus cd found in Fig. 2 of Ref. 331] and Fig. 2 of Ref. 49|. In our case the plots are exactly 
linear due to our perturbative treatment of the three-nucleon interactions. In Fig. H] we 
show i?triton(-^) versus L measured in physical units at LO, NLO, and NNLO. For the NNLO 
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FIG. 4: Triton energy versus periodic lattice length at LO, NLO, and NNLO. For the NNLO 
results we show data for cd = 1.0 fitted to the physical triton energy. 



calculation we show data for cd = 1-0 fitted to the physical triton energy. We see that 
the NLO and NNLO corrections each appear small enough to be treated using perturbation 
theory. 



B. Neutron-deuteron scattering 



Liischer's formula relates the energy levels for a two-body system in a finite periodic 



cube to scattering phase shifts at infinite volume [50|, |51i- We use this method to calculate 
neutron-deuteron scattering phase shifts in the spin-doublet and spin-quartet channels. For 
cubic lattice lengths L = 4,5, 6, 7, 8 we measure the three-nucleon energy levels relative to 
the threshold energy for a non-interacting neutron and deuteron in the same volume. 

In Fig. [5] we plot p cot 6 versus in the center-of-mass frame for the spin-doublet channel. 
Using the pseudo -Etriton constraint for ce, we show NNLO results for cn = —6.0, 0.0, 6.0. 



This can be compared with the physical scattering length 



0.007th. fm 



nd 



-0.645 ± 0.003cxp. ± 



52l | and pseudo scattering length "^and = —0.45 



4) fm resulting from adjusting 



the strength of nn interactions to match np interactions 33|]. A detailed discussion of 
isospin-breaking contributions to the neutron-deuteron scattering lengths can be found in 
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FIG. 5: NNLO results for spin-doublet neutron-deuteron scattering. We plot ;5cot5 versus in 
the center-of-mass frame. 



Ref. 531]. We match to the pseudo scattering length and find that cd hes in the range from 
—6.0 to +6.0. This is a rather loose constraint since we already expect cd ~ 0{1) based on 
the natural size of coefficients under renormalization group transformations. We consider 
alternative methods for constraining cd later in our discussion. 

In Fig. E] we plot pcot6 for the spin-doublet channel for a wider range of p"^. For the 
NNLO calculation we show data for cd = 1-0 fitted to the physical value for -E'triton- The 



experimental results are nd and pd scattering data from the partial wave analysis in Ref. 54 1- 
The dashed line shows an empirical model introduced in Ref. 5J] with a pole singularity 
in pcot6 just below zero energy. This empirical model also accommodates data points for 
the triton and '^He bound states at negative p^. In our lattice data we also find non-trivial 
scattering behavior just below zero energy. The interpretation of these results and possible 
connections with the Efimov effect at finite volume are currently being studied jssl. l56l. Is?]. 
The deuteron break-up threshold is near p"^ = 0.07 fm~^, and the agreement between lattice 
and experimental results for nd scattering is quite good below break-up. Above the break- 
up threshold our analysis using Liischer's finite volume formula does not take into account 
mixing between nd and nnp three-nucleon states. Therefore we expect significant errors in 
the case of strong mixing. 
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FIG. 6: Plot of pcotS versus for the spin-doublet channel in the center-of-mass frame. For 
the NNLO calculation we take cd = 1-0 fitted to the physical value for -Etriton- The experimental 
results are from Ref. 54]. 



In Fig. [7| we plot p cot S for the spin-quartet channel. Again we show experimental results 



for nd and pd scattering from Ref. 
c/) = 1.0 fitted to the physical value for -Etriton 



54j . and for the NNLO calculation we present data for 
The lattice data agree with experimental 



results for nd scattering below deuteron break-up. However significant deviations appear 
above the break-up threshold. This may indicate mixing effects between nd and nnp states. 

V. TRANSFER MATRICES WITH AUXILIARY FIELDS 



For systems with more than three nucleons, sparse-matrix calculations using the lattice 
transfer matrix are not practical at large volume. Instead we use projection Monte Carlo 
with auxiliary fields. A review of the auxiliary-field formalism can be found in Ref. 46] . 
We define M^'^*\7i'j, s, sj) as the leading-order auxiliary-field transfer matrix at time step nt, 



M("*)(7r;,s,S7) =: exp <^ -Hf,,,at 



n,S,I 



nJ 
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FIG. 7: Plot of pcot(5 versus p"^ for the spin-quartet channel in the center-of-mass frame. For 
the NNLO calculation we take cd = 1-0 fitted to the physical value for £^triton- The experimental 
results are from Ref. [54]. 



We can write Mlo as the normalized integral 



J Dn'jDsDsi e 



UtTjUsiJsj e 



where S"!"'^ is the piece of the instantaneous pion action at time step n^, 



n,I 



and S'lr*'' is the auxiliary-field action at time step n^, 



n,n' 



I n,n' 



with 



(59) 



(60) 



St^ = ^ E ^(^' ^t)f-\n - n')sin', nt) + lj2Yl ^^(^' ^t)f-\n - n')siin\ n,), (61) 



(62) 
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The NLO and NNLO interactions are treated using perturbation theory. We let 

n n,S n,S 

n,S,S' n,l 

n,l,S n,l n,l,S 

(63) 

and 

n,I n,S,I n,S,I 

n,S,S',I n,lJ 

n,l,S,I n,l,I n,l,S,I 

(64) 

With these extra fields and linear functionals we define 
= : exp ^ -Hi,^^at - ^^"L ^ Agvr j (n, nt)psf (n) 



nj 



(65) 



In Ref. m 



46 1 there were factors of ^/a't multiplying f/("*)(e) and ?/]"*'*(£). We have 
removed these factors here as they complicate the discussion of the three-body interactions 
at NNLO. Let M^'^^\e) be the normalized integral, 



j Dn'DsDsi e-^^"''-^-*'M("')(7r;, s, si, 



M("*)(£) = ^ . (66) 

/ Div'DsDsj e~'^'i''^*'~'^^^'' 

When all e fields are set to zero we recover Mlo, 

M("*)(0) = Mlo- (67) 
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To first order in perturbation theory the NLO interactions in Mnlo can be written as a 
sum of bihnear derivatives of M^'^*\e) with respect to the e fields at £ = 0, 



S 



3^ 5ep{n,nt) 5ep{n,nt) 



£=0 



2 ^ Ssp(n, rit) Ss^2p(n, rit) 
Similarly for NNLO we have 



£=0 



(68) 



-Mnnlo = -Mnlo + 



where 



(69) 



r(3JV)K) 



(3iV){nO 



(70) 



The three-nucleon contact interaction is 



M 



(mint) / / 



contact 



tt') = -^-DcontactO;* ^ 



6 



1 3 



fcp(n, nt) 



M("*)(£) 



£=0 



and the one-pion exchange interaction has the form 



(3iV)(nt) 



tt') = -DopE-^ ^ A57rj(n,nt) 



OPE 



n,S,I 



5£ps j{n, nt) Sep{n, rit) 



M("')(£) 



£=0 



The three two-pion exchange terms are 



Mr 



(3iV){nt) 



at 



TPEl K^') — ~-DtPE1 — 



[As7r^(n,nt)As7rj(n,nt) - {As7i'j{n,nt)As7i'i{n,nt))] 



n,S,I 



Sep{n, nt) 



£=0 



(71) 



(72) 



(73) 



(3JV)K) 



TPE2 



TPE2" 



X 



p7r/(n, nt)n7ry(n,nt) - (^^(n, nt)n7rj(n, nt))] 



5M("')(£) 
fep(n, nt) 



£=0 



(74) 
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^TPE3 ) = --DtPE3 — 



^ J2 ^s^,S2,S3£h,h,i-Asyh{n,nt)Asyj^in,nt) 



Sep {n, nt) 



■ (75) 

e=0 



n,Si,S2,S3 Il,l2,h P^3-l3 

We extract the properties of the ground state using Euchdean-time projection. Let 
l^freeN^ bc a Slater determinant of free-particle standing waves on the lattice. We construct 
the trial state |^E'(t')) using 

|vi/(t')) = (Msu(4v)^'°|^^^^^>, (76) 

where t' = Lt^at and Lf^ is the number of "outer" time steps. As the notation suggests, the 
transfer matrix Msu(4)/ is invariant under an exact Wigner SU(4) symmetry and acts as an 
approximate low-energy filter. The amplitude Z{t) is defined as 

Z(t) = (vI/(t')|(MLo)^*M*(0), (77) 
where t = Lt-at and Lt- is the number of "inner" time steps. The transient energy 

Ei^o{t + at/2) (78) 
is given by the ratio of the amplitudes for t and t + a*, 

(79) 



p-i?Lo(t+«t/2)-at ^ ^(^ + ^t) 

z{t) ■ 

The ground state energy -Eq.lo equals the asymptotic limit, 



Eo,Lo = lim^Lo(t + «t/2). (80) 



We calculate these Euclidean-time projection amplitudes using auxiliary fields. For a given 
configuration of auxiliary and pion fields, the contribution to the amplitude Z{t) is propor- 
tional to the determinant of an A x A matrix of one-body amplitudes where A is the number 
of nucleons. Integrations over auxiliary and pion field configurations are computed using 



hybrid Monte Carlo. Details of the method can be found in Ref. [4, |46|, |58|, |59 |. 

For the ground state energy at NLO and NNLO we compute expectation values of Mlq; 
Mnlo, ^nnlo inserted in the middle of a string of Mlo transfer matrices, 

ZM,o(i) = (^(^')l (Mlo)^*^/' Mlo {M^^of^^'^ |^(t')) , (81) 
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ZM.,oit) = (Mlo)^'^/' Mnlo (Mlo)^*'/' l^(t')) , (82) 

^M,NLo(^) = (Mlo)^*^/' Mnnlo (Mlo)""'/' l^(t')) • (83) 

From the ratio of amplitudes, 

^^^^ = 1 - AE^^o{t)a, + ■ ■ ■ , (84) 

we define the transient NLO energy correction Aii^NLo(^)- The elhpsis denotes terms which 
are beyond first order in the NLO coefficients. The NLO ground state energy -Ecnlo is 
calculated using 

^o,NLO = ^o,LO + hm AENLo(t)- (85) 

Similarly at NNLO we have 



2'a/lo (t) 

and 



AENNLo(t)at + ■■■ , (86) 



-£^0,NNLO = -Eo,NLO + 1™ AEnNLo(^)- (87) 



VI. PRECISION TESTS OF MONTE CARLO SIMULATIONS 

We use the three-nucleon system to test the auxiliary-field Monte Carlo simulations. The 
same observables are calculated using both auxiliary-field Monte Carlo and the exact transfer 
matrix without auxiliary fields. We choose a small system so that stochastic errors are small 
enough to expose disagreement at the 0.1% — 1% level. We choose the spatial length of the 
lattice to be L = 3 and set the outer time steps Lt^ = and inner time steps Lt- = 4. With 
2048 processors we generate a total of about 10^ hybrid Monte Carlo trajectories. Each 
processor runs completely independent trajectories, and we compute averages and stochastic 
errors by comparing the results of all processors. 

We choose to be a spin-doublet isospin-doublet state built from the Slater deter- 

minant of standing waves IV'2), iV^s) with 

(0| aij(n) \^pi) oc SifiSj^i, (0| aij{n) |?/;2) oc 1, (0| aij{n) {ips) oc 5i,o5i,o- (88) 

In Table [IT] we show Monte Carlo results for the energy (MC) versus exact transfer matrix 
calculations (Exact) at LO, NLO, and NNLO. The NNLO data uses cd = 1.0 with ce 
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TABLE II: Monte Carlo results for the energy versus exact transfer matrix calculations 



Energies 


MC 


Exact 


^Lo(t + at/2) [MeV] 


-48.873(18) 


-48.8823 


A^NLoCt) [MeV] 


0.5509(8) 


0.55100 


A^NNLo(t) [MeV] 


-0.967(3) 


-0.96718 



TABLE III: Monte Carlo results versus exact transfer matrix calculations for the derivative of the 
energy with respect to NLO coefficients 



NLO energy derivatives 


MC 


Exact 


a[AC) MeV'^] 


Q nnQ'7/' 1 o\ 
O.y03/(lzj 


3.9UZZD 




-4.847(2) 


-4.84331 


d(AE^^o(t)) [109 MeV^] 


-2.0105(6) 


-2.01059 


9(^^^^^oit)) MeV^] 


2.9230(14) 


2.92424 


9(f;^^^oit)) Mgg MeV^] 


0.1860(12) 


0.18411 




5.094(2) 


5.09371 


d{AE^^o,{t)) J1Q9 MeV^] 


-1.5892(3) 


-1.58898 


diAE^^oit)) [109 MeV5] 


6.8019(11) 


6.80197 


diAE^^^om mq9 MeV^] 


0.3417(2) 


0.34164 


9{Ai?HLoW) [109 MeV5] 


-1.0092(5) 


-1.00932 



fitted to the physical triton energy. In Table IIIII we compare Monte Carlo results with 
exact transfer matrix calculations for the derivative of the energy with respect to each NLO 
coefficient. In Table IIVI we make the same comparison for the derivative of the energy 
with respect to each NNLO coefficient. The numbers in parentheses are the estimated 
stochastic errors. In all cases the agreement between Monte Carlo results and exact transfer 
calculations is consistent with estimated stochastic errors. 
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TABLE IV: Monte Carlo results versus exact transfer matrix calculations for the derivative of the 
energy with respect to NNLO coefficients 



NNLO energy derivatives 


MC 


Exact 


C(, -(^contact ) ^ ' 


1.162(4) 


1.1609 


a(AW^ [10^ MeV5] 


-5.858(6) 


-5.8623 


9(AWoM) mqS MeV^l 

O(U'YPEl) L J 


14.46(5) 


14.468 


5(ABnnloW) mqS MeV^l 
o(Dtpe2) l j 


2.24(3) 


2.2595 


5(ABnnloW) mqS MeV^l 

o(Dtpe3) ^ J 


-10.02(6) 


-10.022 



VII. ENERGY OF THE a PARTICLE 

We simulate the a particle on cubic periodic lattices with length L = 5, 6, 7, 8. These 
correspond with physical lengths L = 9.9, 11.8, 13.8, 15.8 fm. For jvl/^'''^''^ we take the Slater 
determinant formed by standing waves 

(0| aij{n) oc 5ifl5j^i, (0| aij(n) \tp2) cx (89) 

(0| aij(n) 1^3) oc Si^oSjfl, (0| aij{n) \tp^) oc (5i,i5j,o- (90) 

For each value of Lf. a total of about 6 x 10^ hybrid Monte Carlo trajectories are generated by 
2048 processors, each running completely independent trajectories. Averages and stochastic 
errors are computed by comparing the results of all processors. 
For L = 5, 6, 7, 8 we show Monte Carlo results in Fig. [8] for 

^LOW, ^LOW+A^NLOW, ^LOW+A^NLOW+A^NNLOW, (91) 

versus Euclidean time t. These are labelled as LO, NLO, and NNLO respectively. The 
NNLO data uses cd = 1-0 with ce fitted to the physical triton energy. In addition to the 
Monte Carlo data we plot the asymptotic forms for each term in Eq. (19T!) using 

ELo(t) ~ Eo,LO + ^Loe-'^-*, (92) 

AENLo(t) ~ ^0,NLO - Eo,LO + ^NLoe"''' */', (93) 
A£JnNLo(^) ~ -E'o.NNLO — -Eo.NLO + -SnNLOC"'^^'*''^. (94) 

23 



L = 5 



L = 6 




FIG. 8: Monte Carlo results for a-particle energies versus Euclidean time t at LO, NLO, and 
NNLO. We also plot fitted asymptotic expressions. The NNLO data uses cd = 1-0 with ce fitted 
to the physical triton energy. 



The unknown coefficients Alo, -Snlo; -BnnlO) and energy gap 5E, are determined by least 
squares fitting. The e~^^'^ dependence in Eq. fl92l) comes from the contribution of low-energy 
excitations with energy gap SE above the a particle. The e"*^^ */^ dependence in Eq. ( p3|) 
is due to matrix elements of AfNLO between the a particle and low-energy excitations at 
5E. Similarly the e"''-^'*/^ dependence in Eq. (1941) is from matrix elements of Mnnlo- The 
chi-squared per degree of freedom for the fits are 1.5 for L = 5, 1.0 for L = 6, 0.7 for L = 7, 
and 0.6 for L = 8. The small relative size of the NLO and NNLO energy corrections suggest 
that the perturbative treatment of these terms appears reasonable. 
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FIG. 9: Energy of the a particle at NNLO versus c/j. The contact interaction ce is fitted to the 
physical triton energy. The dotted line is the estimated Coulomb-subtracted energy —29.0 MeV. 



In Fig. Owe plot the NNLO a-particle energy versus c^, with ce fitted to the physical 
triton energy. The bands indicate the estimated error due to stochastic noise and asymptotic 
fits at large t. The a energy should approach the infinite volume limit from below, similar 
to our results for the triton energy. Hence the deviation between data at L = 7 and L = 8 is 
likely due to stochastic noise and fit errors rather than finite volume effects. The a energy 
shown at —29.0 MeV is the estimated Coulomb-subtracted energy jsO]. At large volumes the 
best agreement with the Coulomb-subtracted a energy occurs at ~ —4. The a binding 
increases in strength by 0.2 MeV for each unit increase in cd, and so we find reasonable 
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agreement for all values of cd ~ 0(1). If instead we fit c^; according to the pseudo triton 
energy, then the lines in Fig. [9] shift downward in energy by about 2 MeV. The pseudo a 
energy with pp- and nn-forces matched to np-forces is estimated to be —29.8(1) MeV |33l ]. 



VIII. SUMMARY AND DISCUSSION 

We have presented the first study of low-energy few-nucleon systems on the lattice at 
next-to-next-to-leading order in chiral effective field theory. We computed nucleon-nucleon 
phase shifts and S-D mixing angle on the lattice and used scattering data to determine 
unknown two-nucleon operator coefficients. In the three-nucleon system we calculated the 
triton energy and determined neutron-deuteron phase shifts using Liischer's finite volume 
method. These were used to constrain the two cutoff-dependent three-body coefficients, 
and ce- For the four-nucleon system we recast the lattice action in terms of auxiliary fields 
and used projection Monte Carlo to calculate the energy of the a particle. 

There are several ways in which the lattice calculations presented here can be improved 
and extended in future work. One improvement is the inclusion of isospin-breaking ef- 
fects due to Coulomb interactions and quark mass differences. The framework for isospin- 
violating effects in chiral effective field theory has been developed over the past decade 



60, 



61 



62 



63 



64 



65 



67 



681]. Work is currently underway to implement these effects 
within the lattice formalism. 

Another area of improvement concerns the P-wave phase shifts for our leading-order lat- 
tice action. NLO corrections to the phase shifts are substantial for nucleon momenta above 
100 MeV. This seems not to cause any problems for the light 5'-shell nuclei considered here. 
However for P-shell nuclei we may find corrections strong enough to spoil the perturbative 
treatment of higher-order effects. In Ref. [3] this problem has already been resolved in 
neutron matter calculations using a new leading-order action LO3. This lattice action uses 
spin-isospin projection operators to produce Gaussian smearing only in even partial wave 
channels. The implementation of the LO3 action in Monte Carlo simulations with both 
protons and neutrons is computationally more intensive than the pure neutron simulations 
in Ref. The auxiliary-field formalism requires a total of 16 auxiliary fields and some 
increase in the sign/phase cancellations relative to the LO2 action. Nevertheless studies of 
light P-shell nuclei using the LO3 action are planned in the near future. 
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One recent paper constrains the cutoff-dependent coefficient cd from the triton beta decay 
rate \^^. From the point of view of computing the spectrum of hght nuclei, the easiest 
method for determining cd is by means of the a-particle energy. This has the added benefit 
of removing systematic errors from the four nucleon system. If however we also want to 
accurately describe the chiral interactions of nucleons and light nuclei with soft pions, then 
it would be best to measure cd directly from the SU(2) axial coupling to two nucleon states. 



This is studied in Ref. 



701 ] using pion production data in pp scattering. Unfortunately 



the pion production energy threshold is too high to be accessible at our lattice spacing and 
extrapolations are required. In the future another approach may be possible using direct 
theory-to-theory matching. In this technique lattice QCD would be used to calculate the 
axial charge of two-neutron scattering states in a periodic cube. This calculated value of 
the axial charge could then be used to fix cd for any chosen lattice spacing in lattice chiral 
effective field theory. 
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APPENDIX: LATTICE NOTATION 



The vector n represents integer-valued lattice vectors on a three-dimensional spatial lat- 
tice, and p, q, k represent integer- valued momentum lattice vectors. Z = 1, 2, 3 are unit 
lattice vectors in the spatial directions, a is the spatial lattice spacing, and L is the length 
of the cubic spatial lattice in each direction. The lattice time step is at, and rit labels 
the number of time steps. We define at as the ratio between lattice spacings, at = at/ a. 
Throughout our lattice discussion we use dimensionless parameters and operators, which 
correspond with physical values multiplied by the appropriate power of a. Final results are 
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presented in physical units with the corresponding unit stated exphcitly. 

We use a and to denote annihilation and creation operators. We make explicit all 
spin and isospin indices, 

O'lfi = 0'l,p: 0,1,1 = (^l,n- 

The first subscript is for spin and the second subscript is for isospin. We use tj with 
/ = 1,2,3 to represent Pauli matrices acting in isospin space and as with S = 1,2,3 to 
represent Pauli matrices acting in spin space. 

We use the eight vertices of a unit cube on the lattice to define spatial derivatives. For 
each spatial direction / = 1, 2, 3 and any lattice function f{n), let 

^if(^) = \ E i-'^r^'fin + u), z7 =z/ii + z/22 + z/33. (A.3) 

We also define the double spatial derivative along direction /, 

/(n) = f{n + + f{n - 1) - 2/(n). (A.4) 
For the three-body NNLO interactions we also use the notation 

°f(^) = l Yl /(^+'^)' z7 =z/ii + z/22 + z/33. (A.5) 

1. Local densities and currents 

We define the local density, 

P'^'^in) = Yl 4,,in)a^An), (A.6) 

i,j=0,l 



which is invariant under Wigner's SU(4) symmetry 39j. Similarly we define the local spin 
density for = 1, 2, 3, 

Pl'''i^)= 4,j(^)\-^s]ii>ai'j{n), (A.7) 

i,j,i'=0,l 

isospin density for / = 1,2,3, 

pf'''{n)= «l(^)N,ya..'(^), (A.8) 

i,j,j'=0,l 
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and spin-isospin density for 5", / = 1, 2, 3, 

Psfi^)= Yl ^h(^)l^s]iAri]jj,ai,j,{n). (A.9) 

i,j,i' ,j'=0,l 

For each static density we also have an associated current density. Similar to the defini- 
tion of the lattice derivative A; in Eq. flA.Sp . we use the eight vertices of a unit cube, 

i7=z/ii + z/22 + z/33, (A. 10) 

for Ui, 1/2, 1^3 = 0, 1. Let for / = 1, 2, 3 be the result of reflecting the /*'^-component of 

V about the center of the cube, 

= z7+ (1 - 2z/i)f. (A.ll) 

Omitting factors of i and 1/m, we can write the /*'^-component of the SU(4)-invariant current 
density as 

nf '"(n) = \ Y. E + (A.12) 

Similarly the /*'^-component of spin current density is 

</(^) = i E E {-ir^'al{n + V{-l))[asl,a,,{n + iJ), (A.13) 

'^l,l'2,'^3=0,l i,j,j'=0,l 

/*^-component of isospin current density is 

y\,u2,v3=0,l i,j,j'=0,l 

and /^'^-component of spin-isospin current density is 

nfc(^) = ^ E E (-l)'^'+H,(^ + ^(-0)k5]..Kya.'.'(^ + ^)- (A.15) 

2^1,1^2, 1^3=0,1 ,j'=0,l 

2. Instantaneous free pion action 

The lattice action for free pions with purely instantaneous propagation is 
S'ttttIvt/) = + 3) ^ 7r/(n,nt)7r/(n,ni) - Of ^ 7ii{n,nt)'Ki{n + l,nt), (A.16) 

n,nt,I n,nt,I,l 

where tt/ is the pion field labelled with isospin index J. It is convenient to define a rescaled 
pion field, vr^, 

7i'j{n,nt) = ^/q^7^■I{n,nt), (A.17) 
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= at{ml + 6). 

Then 

In momentum space the action is 



n,nt,I,l 



I,k 



1 Q;^ 

2 " ^ 



cos ki 



(A.18) 



(A.19) 



(A.20) 



The instantaneous pion correlation function at spatial separation n is 
//^ X //r; A /i:'7r;7r;(n,nt)7r;(0,nt) exp[-5^^] 



J Dn'j cxp [—St, 



-ik-n 



(A.21) 



where 



DM 



1 - — E; COS k 



It is also useful to define the two-derivative pion correlator, Gs^Sjl^)? 



(A.22) 



Gs^S2 (^) = i^syiifi, nt)As27rj{Q, Ut) ) (no sum on /) 



1 
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(A.23) 
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